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FORTRAN PROGRAM FOR GENERATING A TWO-DIMENSIONAL ORTHOGONAL 
MESH BETWEEN TWO ARBITRARY BOUNDARIES 
by William D. McNally 
Lewis Research Center 

SUMMARY 

A FORTRAN IV program is presented which computes and plots the coordinates for 
a two-dimensional orthogonal mesh in the region between the walls of a flow channel. 

The program is designed for a channel containing a body about which flow passes and 
which spans the channel from one wall to the other. However, the condition that the 
channel contain an immersed body can be easily removed from the program. 

Input to the program is brief. It consists of arrays of point coordinates describing 
the two channel walls and the front and rear edges of the submerged body. Information 
is also given for the number of orthogonal mesh lines desired in the throughflow and 
normal directions. Output consists of the coordinates of the generated orthogonal mesh 
points as well as the angles of the mesh at these points with the horizontal plane. One 
of the routines in the program also generates a microfilm plot of the channel and the 
submerged body, as well as a plot of the orthogonal mesh. All routines in the program 
except the plot routine are in general FORTRAN IV code and could be easily transferred 
to other computing equipment. The plot routine makes use of Lewis in-house subroutines 
and equipment and would require recoding at another facility. 

The geometry of each of the channel walls input to the program must be smooth 
enough that it can be easily fit with a cubic spline curve; that is, it should not contain 
any sharp discontinuities in slope. The space between these walls is divided into equal 
increments across the channel, and spline curves through these points constitute the 
M horizontal M , or throughflow, orthogonal lines. The ’’vertical", or normal, orthogo- 
nals are obtained by means of a simple "predictor- corrector" type process in which 
normals to adjacent orthogonal lines are used. The program is restricted to channels 
which are essentially horizontal; that is, at no point should either of the channel walls 
be inclined more than 45° to the horizontal. 

This report includes a listing of the program, with an explanation of the input re- 
quired and the output generated. Numerical examples are also included. Running times 
are under 0. 1 minute on IBM 7094 or 360-67 equipment. 
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INTRODUCTION 


Finite difference solutions to partial differential fluid flow equations require that a 
finite difference mesh be placed upon the solution region. And in most cases, the terms 
in the finite difference equations will be simplified if this mesh is orthogonal. This 
report describes a FORTRAN IV program which computes and plots coordinates for a 
two-dimensional orthogonal mesh in the region between the walls of a flow channel. 

The primary application for this program at NASA Lewis Research Center is for 
analytical solutions for flow in turbomachinery passages containing rotor or stator blade 
rows. For this reason, the program is designed for a channel containing a solid body 
about which the flow passes and which spans the channel from one wall to the other. The 
input to the program is designed so that the mesh spacing can be varied in the regions 
upstream of the immersed blade row, in the general space occupied by the blade row, 
and downstream of it. In this way a region of particular interest can be covered with a 
finer mesh than is used in the other two regions. 

This "immersed body" feature of the program can be easily removed, however, 
since the only way the body geometry affects the generated mesh is by designating the 
points of separation of the three mesh-size regions. The resulting program would have 
wide application in many fields where finite difference solutions are desired between the 
walls of an arbitrarily shaped channel . 

Input to the program consists primarily of spline points describing the channel walls 
and the edges of the submerged body. These arrays of points must be smooth enough 
that they are capable of being accurately fitted with cubic spline curves . Therefore, 
they cannot contain any sharp discontinuities in slope. 

The program divides the space between the walls into equal increments across the 
channel. Spline curves are then fit through the resulting points to obtain the throughflow, 
or streamwise or "horizontal", orthogonals. The normal, or "vertical", orthogonals 
are obtained by means of a simple "predictor- corrector" technique in which normals to 
adjacent orthogonal lines are used. This technique is analogous to the second- order 
Runge-Kutta method for solving ordinary differential equations known as the improved 
Euler method or Heun's method (ref. 1). 

The program is restricted to channels which are essentially horizontal; that is, the 
channel walls should not be inclined by more than about 45° with the horizontal plane. 

For axial turbomachinery this is not a serious restriction. This limitation could be 
removed in a more general recoding of the program. The program could also be ex- 
tended to two-dimensional geometries with walls having noncontinuous slopes, as well 
as to some three-dimensional geometries. 

Output consists of printed coordinates of the generated orthogonal mesh points, as 
well as the angles of the mesh lines at these points with the horizontal plane. One of the 
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program's subroutines generates a microfilm plot of the channel input geometry and of 
the calculated orthogonal mesh. All routines in the program except the plot routine are 
in general FORTRAN IV code and could be easily transferred to other computing equip- 
ment. The plot routine makes use of Lewis in-house subroutines and equipment and 
would require recoding at another facility. Run times for average mesh sizes are under 
0. 1 minute on IBM 7094 or 360-67 equipment. 

This report includes a description of the method used to generate the orthogonal 
mesh, as well as a description of each of the subroutines in the program. A complete 
listing is also included. The input and output are thoroughly described, and numerical 
examples are given. The Lewis plot routines are described in the appendix. 


DESCRIPTION OF METHOD 

The program makes use of four principal subroutines - INPUT, MESH, OUTPUT, 
and OUPLOT. The bulk of the work is done in MESH, which uses four auxiliary 
routines - ROOT, CROSS, SPLINE, AND SPLINT. About 400 lines of code are used in 
these eight routines. The calling arrangement is illustrated in figure 1. 



Figure 1. - Calling sequence of subroutines. 


The INPUT routine reads in data, the OUTPUT routine prints results, and OUTPLOT 
makes a microfilm plot of the generated mesh. So the bulk of the mesh generation 
algorithm is contained in MESH, which is described in the following paragraphs. 

The mesh generation algorithm begins with input data similar to those shown in 
figure 2. Two sets of coordinates (ZBOT,RBOT and ZTOP,RTOP) describe the upper 
and lower channel walls, and two others (ZLE,RLE and ZTE,RTE) describe the front 
and rear of the confined body. The coordinates defining the body are not used as part 
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Figure 2. - Schematic representation of input variables. 


of the mesh generation algorithm. It is the shape of the channel walls which fully de- 
termines the orthogonal mesh. The intersections of the edges of the body with the 
bottom channel wall (ZLE(l)and ZTE(l)) are the points which separate the three regions 
of varying mesh size. These regions correspond generally to the areas upstream of the 
body, on the body, and downstream of it. The Z- coordinates, ZBEGIN and ZEND, de- 
fine the extremes of the orthogonal mesh where its left and right boundaries intersect 
the bottom channel wall. These two parameters, along with the number of mesh spaces 
requested in each region, determine the mesh spacing upstream and downstream of the 
central body. The two points where the left and right mesh boundaries intersect the 
upper channel wall must remain arbitrary, since they will be defined by the mesh gener- 
ation algorithm. 

The MESH subroutine begins by extending lines vertically from each of the input 
points on the lower boundary to the upper channel wall (see fig. 3). Each of these lines 
is then divided into equal increments, the number depending upon the number of 
mesh lines desired in the horizontal, or streamwise, direction. The resulting array of 
points along the radial lines is called RRAD. Spline curves are fit through the points in 
the RRAD array to give the horizontal orthogonals. These are shown in figure 4. 

The procedure for calculating the vertical orthogonal links between adjacent hori- 
zontal orthogonals is illustrated in figure 5. This procedure is analogous to the tech- 
nique for solving ordinary differential equations known as the improved Euler method 
or Heun’s method (ref. 1). In this technique, the solution at the next unknown point is 


4 


I 






first predicted from the solution at a known point by using the Euler method. A second 
solution is then calculated by using a new slope obtained from the predicted point. The 
final solution is a corrected combination of these two solutions. 

In this program, a similar process is used to construct the vertical links between 
adjacent horizontal orthogonals. These links are constructed one at a time, moving 
from left to right between each pair of adjacent horizontal orthogonals (see fig. 6). 

When all links are calculated for one pair of orthogonals, the algorithm moves vertically 
to the next pair and again constructs links, moving from left to right. 

Each link begins at a known (previously calculated) mesh point on the lower hori- 
zontal orthogonal with coordinates REFZ, REFR. A normal to this point is constructed 
(line (T) in fig. 5), and the intersection (ZNL,RNL) is calculated on the upper ortho- 
gonal curve. The subroutines ROOT and CROSS are used to locate this intersection. 
(This process will be explained in more detail in the following paragraph. ) The slope 
of the upper orthogonal at the point of intersection (ZNL,RNL) is computed; this is 
SLU. Line ( 2 ) in figure 5 is then constructed in such a way that it is perpendicular to 
the SLU line and also passes through the reference point (REFZ, REFR). The coordi- 
nates of the intersection of line ( 2 ) with the SLU line are calculated; these are ZNU 
and RNU. This point is a small distance away from the upper orthogonal curve due to 
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Figure 6. - Process for generating "vertical" orthogonal links. 

the change in slope of that curve in moving from point ZNL,RNL to point ZNU, RNU. 

The final location of the new orthogonal mesh point on the upper curve has a 
Z-coordinate (ZOM) which is the average of ZNL and ZNU. The corresponding 
R- coordinate (ROM) is calculated by spline interpolation, using the value of ZOM and 
coordinates in the RRAD array. The process of constructing vertical links is continued 
until the upper channel wall is reached by all vertical orthogonals. This completes the 
generation of the mesh. 

Let us return to the procedure by which the intersection of line (T) in figure 5 with 
an upper horizontal orthogonal curve is calculated. The MESH subroutine establishes 
Z-coordinate bounds (BNDL and BNDR) on the upper orthogonal curve between which the 
intersection is certain to occur (see fig. 7). These bounds are calculated by using the 
slope of line (T) and an approximate distance between the two horizontal orthogonals. 

One of these bounds is always equal to REFZ (see fig. 7). The other is established by 
estimating how far the intersection with the upper orthogonal is away (in the Z-direction) 
from REFZ and then doubling that distance for a factor of safety. Once the two bounds 
are set, ROOT is called to find the intersection. 

The ROOT subroutine calculates the intersection of a straight line with a curve 
(the upper horizontal orthogonal in this application), it does this by systematically 
choosing Z- coordinates, calculating the distance in the R-direction between the straight 
line (line (T) in fig. 7) and the orthogonal curve at each of these coordinates, and con- 
verging to the intersection. ROOT begins with a Z-coordinate halfway between the bounds 
BNDL and BNDR and computes the distance between the corresponding R- coordinates of 
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line @ and the upper curve at the given value of Z. This distance is called RMR 
(see fig. 7). Based on the sign of RMR, ROOT decides on which side of the chosen 
Z-coordinate the intersection lies, and moves the left or right bound inward accordingly. 
The resulting interval is again cut in half and a new RMR computed. This process of 
dividing the remaining interval is repeated 20 times by the ROOT subroutine. Each 
time the remaining interval is cut in half, so that after 20 iterations the intersection is 
accurately located. 

CROSS is the routine called 20 times during each call on ROOT. CROSS calculates 
RMR for a given value of Z (see fig. 8). On the straight line (line (7) in fig. 8) normal 
to the lower orthogonal, the R- coordinate (RLINE) is computed for the given value of Z. 
For the same Z, the R- coordinate (RCURV) and slope (SL) are computed on the upper 
orthogonal curve. Note that RLINE can be of greater magnitude than RCURV, as it is 
for the dashed-line value of Z in figure 8. RMR is computed as the difference between 
RCURV and RLINE. 

An example of a completed orthogonal mesh is shown in figure 9. The body geome- 
try has once again been superimposed on the mesh region in this figure. Notice that it 
did not affect the generation of the orthogonal mesh, except that the mesh spacing has 
been varied in the regions forward of the body, on the body, and to the rear of the body. 
This spacing is controlled by the input. 


8 




USE OF PROGRAM 


The orthogonal mesh program is very easy to use. Input consists of a small number 
of coordinate arrays describing boundary geometry. These arrays are very easy to 
assemble. Following sections describe the input and output thoroughly. Numerical ex- 
amples are given to illustrate what the program does and to give an overview of the in- 
put. 


Numerical Examples 

Six examples are given which illustrate the range of geometries for which meshes 
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Z DIRECTION 


ORTHOGONAL MESH 


Figure 10. - Converging channel with fine mesh. 


can be generated. The first of these examples is the converging channel which has been 
illustrated in previous figures in the section DESCRIPTION OF METHOD. This channel 
was shown in figure 9 with a coarse mesh. It is illustrated in figure 10 with a more 
commonly desired fine mesh suitable for finite difference calculations. The complete 
input for this example is listed in table I, and the output is shown and discussed in the 
section Output. 
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TABLE I. - INPUT FOR CONVERGING- CHANNEL EXAMPLE 


ORTHOGONAL 

MFSH 

TEST 

CASE- FINE 

MESH 




MESH INPUT 

DATA 

ANC 

INPUT ARRAY 

BOUNDS 




ML 

MT 

MZ 

MR 

NBOT 

NTOP 

NLE 

NT E 

a 

20 

28 

20 

11 

7 

6 

4 

Z BEGIN 



ZEND 






-0.250000 0.750000 


BOUNDARY INPUT DATA 


ZBOT ARRAY 








-0.500000 

-0.300000 

-0.100000 

0 

0.100000 

0.200000 

0.3O0000 

0.400000 

0.500000 

0.700000 

0.900000 






ROOT ARRAY 
0.878000 

0.920000 

0.968000 

1.000000 

1.04^000 

1.1000C0 

1.152000 

1.200000 

1.248000 

1.330000 

1.383000 






ZTOP ARRAY 
-0.500000 

-0. 100000 

0.100000 

0.250000 

0.400000 

0.600000 

0.900000 


RTOP ARRAY 
1.961000 

1.91U000 

1.8B2000 

1.86D000 

1.84000C 

1.825000 

1.808000 


BODY INPUT DATA 








ZLE ARRAY 

0 

0. 230000E-01 

0.390000E-01 

0.470000E-01 

0.50000CE-01 

0 • 500000E-01 



RLE ARRAY 
l .000000 

1.180000 

1.360000 

1.540000 

1.720000 

1.892000 



ZTE ARRAY 
0.400000 

0.346000 

0.317000 

0.300000 





RTE ARRAY 
1.200000 

1.410000 

1.620000 

1.853000 







Z DIRECTION 

BOUNDARIES OP PLOT REGION 



Z DIRECTION 

BOUNDARIES OF PLOT REGION 


Z DIRECTION 

ORTHOGONAL MESH 



ORTHOGONAL MESH 


Z DIRECTION 



Figure 1L - Input and output plots of typical geometries and meshes. 



DIRECTION R DIRECTION 





The microfilm plots of input points and generated meshes for five other channel 
shapes appear in figure 11. 


Input 

Figure 12 shows the placement of input variables on data cards. The first input 
card is for a title, which identifies the data set and is printed on the output. The user 



Figure 12. - Input data format 


may type whatever information he wishes in any of the columns of this card. The re- 
maining cards are for input data. 

Many sets of data may be executed by the program in one run by simply placing them 
in sequence behind the program. 

The input variables are defined in the next section. That section is followed by some 
special instructions for preparing the input data. 
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Dictionary of input v ariables. - The following variables are given as input to the 

program and are illustrated in figure 13. 

ML number of mesh lines between left boundary of orthogonal mesh 

(ZBEGIN) and front, or leading, edge of body (ZLE(l)) 

MR number of mesh lines between bottom and top boundaries of orthogonal 

mesh 

MT number of mesh lines between left boundary of orthogonal mesh 

(ZBEGIN) and rear, or trailing, edge of body (ZTE(l)) 

MZ number of mesh lines between left and right boundaries of orthogonal 

mesh (ZBEGIN to ZEND) 

NBOT number of points in ZBOT and RBOT arrays defining bottom boundary 

of region 

NLE number of points in ZLE and RLE arrays defining front, or leading, 

edge of body 

NTE number of points in ZTE and RTE arrays defining rear, or trailing, 

edge of body 

NTOP number of points in ZTOP and RTOP arrays defining top boundary of 

region 
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RBOT 


array of R- coordinates of input points defining bottom boundary of 
region 

RLE array of R- coordinates of input points defining front, or leading, edge 

of body 

RTE array of R- coordinates of input points defining rear, or trailing, edge 

of body 

RTOP array of R- coordinates of input points defining top boundary of region 

ZBEGIN Z- coordinate of intersection of first vertical mesh line with bottom 

boundary of region 

ZBOT array of Z- coordinates of input points defining bottom boundary of 

region 

ZEND Z- coordinate of intersection of final vertical mesh line with bottom 

boundary of region 

ZLE array of Z- coordinates of input points defining front, or leading, edge 

of body 

ZTE array of Z- coordinates of input points defining rear, or trailing, edge 

of body 

ZTOP array of Z- coordinates of input points defining top boundary of region 

Special instructions for preparing inpu t. - The following instructions and restrictions 

should be observed when preparing input: 

(1) The number of mesh lines in either direction must not exceed 50: 

MZ < 50 
MR < 50 

(2) The number of input coordinates on any surface must not exceed 50: 

NBOT < 50 
NTOP < 50 
NLE < 50 
NTE < 50 
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(3) The origin of coordinates (Z = 0, R = 0) can be placed in any location convenient 
to the user. 

(4) The arrays of coordinates describing body geometry should begin and end with 
points on the channel walls. If this condition is not met, the program will execute; but 
mesh spacing will be altered depending upon where these arrays of coordinates begin and 
end. 

(5) The horizontal orthogonals are generated by means of spline curves through 
points in the RRAD array. The RRAD array, in turn, is a function of the number of 
points given in the ZBOT array of input. If either channel wall has curvature to it, 
enough points should be given in ZBOT (5 to 10 at least) to ensure an RRAD array with 
enough points to guarantee accurate horizontal orthogonals. 

(6) No points on the generated orthogonal mesh should be outside the extremes of 
the input ZBOT and ZTOP arrays. ZBOT and ZTOP should contain values sufficiently 
outside the limits of ZBEGIN and ZEND to prevent this condition from occurring (see 
fig. 14). If the orthogonal mesh does pass beyond the bounds of ZBOT or ZTOP, the 
SPLINE routines will be used for extrapolation; and the ROOT routine may fail to con- 
verge, stopping the program. 

(7) The array ZTOP should contain Z-coordinates which at least cover the same 
range as those in the ZBOT array. Otherwise, the vertical lines extended upward from 
the bottom channel wall (see fig. 3) will have to intersect an extrapolated portion of the 
upper channel wall. The program will execute under these conditions, but the portion of 
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the upper wall generated by spline extrapolation may not conform to what the user would 
input in its place. 


Output 


Output from the program consists of two parts: a computer listing with printed 
tables of orthogonal mesh coordinates and angles, and microfilm plots of both input 
geometry and generated orthogonal mesh. 

A partial listing of the output from the converging channel of the first example is 
given in table n. The microfilm plots from this example are given in figure 10. 


TABLE n. - PARTIAL OUTPUT FOR CONVERGING- CHANNEL EXAMPLE 


************* ********* **** 
*** **> 

*** L CUOKUl NATES UF **> 
♦** ORTHOGONAL MESH **’ 

*** **i 

**************************, 


** BOTTOM ROM NO. 

-o.?5coon 
0.333334E-01 
U.30000U 
c. 618750 

1 ** 

-0.214286 

0.6666676-01 

0.333333 

0.662500 

-0.178571 
0. 1000006 *-00 
0.366667 
0.70b250 

-0.142857 

0.133333 

0.400000 

0.75C0O0 

-3.137143 
C. 156667 
0.443750 

-0.714286E-01 
0. 200000 
0.487500 

-0.357 143E-01 
C. 233333 
C. 531250 

-0.1 86265E-08 
0.266667 
0.575000 

** ROW NO- 2 ** 
-0.260833 
(i> 17641 9F-0 1 
u. 285571 
U. *>09346 

-0.225334 

0.5005616-01 

0.319556 

0.654000 

-0.189948 
0. 8257236-0 l 
0.353388 
0.69 8698 

-0.1546 63 
0.115429 
0.38 7076 
0.743360 

-0.119461 
0.148905 
C. 431257 

-0. 84363 7E-0 1 
0.182907 
0.475628 

-C.494692E-01 

0.217214 

0.520157 

-0.14729IE-CI 

0.251453 

0.564733 

*♦ ROM NO. 3 ** 
-0.270797 
6. 336329b -02 
0.271495 
0.600 201 

-0.235447 
0* 348 1 36E— 01 
0.306189 
0.645746 

-0.200325 

O.664123E-01 

0.340566 

0.691380 

-0.165403 
0. 9 84879E-01 
0.374634 
C. 736963 

-0.130650 

0.131704 

0.419175 

-0. 960822 E-01 
0. 166089 
0. 464085 

-0. 6 19012E-01 
0.201299 
0.509350 

-0. 2 806C8E-01 
0.236525 
0.554738 

** RCIm NU. 4 ** 
-0.279901 
— 0.946350E— 02 
0.257829 
0.591356 

-0.244640 
0.21G250E— 01 
C. 293290 
0.637774 

-0.209716 

0.516590E-C1 

0.328260 

0.684328 

-0.175092 

0.828075E-D1 

0.362737 

0.730831 

-0.140724 

0.115367 

0.407596 

-0.106613 

0.149780 

0.452948 

-0.730324E-01 
0.1 85677 
0.498884 

-f. 399967E-01 
0.221940 
n. 545062 

** ROW NU. 5 ** 
-0.288153 
-0.20B329E-01 
0.244641 
0.582854 

-0.252922 

0.873441E-02 

0.280919 

0.630122 

-0.218135 
C. 384014E-U 
0.316533 
0.677573 

-0.183745 
C. 685735E-D1 
0.351449 
0.724988 

-D. 149696 
0.100214 
C. 396615 

-0.115979 

0.134248 

0.442296 

-0 .8 28917E-01 
0.170515 
0.488822 

-C.505551E-01 

0.207763 

0.535756 

** ROW NU. 6 ** 
-0.295559 
—O. 30 76126— 01 
0.232004 
0.574742 

—0 .260303 
-0.204624t-02 
0.269149 
0.622828 

-0.225595 
0.2668556-01 
G. 305447 
0.671148 

-0.191376 
0.5567696-31 
0.340 836 
0.719459 

-0.157583 
0 • 6650C0E-01 
0.386316 

-0.124195 
C. 119772 
0.432214 

-0. 915096E-01 
0.156050 
C. 479231 

-0.597631E-01 

0.194079 

0.526873 

** ROM NU. 7 ** 
-0.302118 
— 0.392 746F— 0 1 
0.220003 
0.567066 

-0.266787 
-0 • 1 1 3282 E -01 
0.258057 
0.615934 

-0. 232104 
0.165253E— 01 
0.295071 
0.665087 

-0.197999 
0.4476516-31 
C. 330962 
0.714268 

-0.164399 

0.743638E-01 

0.376768 

-0.131275 
0. 106621 
0.422789 

-0.989149E-01 

0.142526 

0.470187 

-0.6765C6E-01 
0.181034 
° • 518469 

** ROM NO. 8 ** 
-0.307831 
—0 . 46403 3F—0 1 
0.208734 
0.559878 

-0.272379 
— 0. 191368E-0 1 
0.247728 
U. 60 9481 

-0. 237668 
0. 791 1 80 1— D 2 
0.285478 
0.659424 

-0.233621 
0. 352528E-D L 
C. 321896 
0.709441 

-0.17C155 
0.63 8 549E— Cl 
0.368036 

-0.137232 
0.950118E-01 
C. 414106 

-0.105131 

C.13C180 

0.461763 

-0.74246BE-01 

0.168830 

0.5106C3 

** ROM NO. 9 ** 
-0.312691 

-0.277080 

-0.242293 

-C. 208250 

-0.174860 

-0.142078 

-0.110172 

-P.795773E-01 
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TABLE n. 


Concluded. PARTIAL OUTPUT FOR CONVERGING- CHANNEL EXAMPLE 


**♦**+***♦***************** 
*•* *** 

*** R COORDINATES OF *** 
*** ORTHOGONAL HESH *** 
*** *** 

*************************** 


** BOTTOM ROM NU. 1 ** 


0.9 3094 8 
1 .0 i 3058 
1.I520Q0 
1.299921 

0.939051 

1.027663 

1.168154 

1.316761 

0.947521 

1.044000 

1.164075 

1.332091 

0.956470 

1.062059 

1.20000D 

1.345625 

0.966011 

1.081053 

1.221124 

0.976274 

1.100000 

1.242108 

0.987504 

1.118113 

1.262418 

1.000000 

1.135385 

1.281749 

** ROW NO. 2 ** 
0.981819 
1.053565 
1.182247 
1.323919 

0.989140 

1,065995 

1.197764 

1,340348 

0.996724 

1.079952 

1.212897 

1.355339 

1.004676 

1.095724 

1.227897 

1.366737 

1.013097 

1.113166 

1.247805 

1.022084 

1.131356 

1.267857 

1.031783 

1.149213 

1.287481 

1.042462 

1.166127 

1.306235 

** ROW NO. 3 ** 
1.033250 
1.095780 
1.713518 
1.348413 

1.039846 

1.106279 

1.228397 

1.364350 

1.046616 

1.116092 

1.24*747 

1.378924 

1.053658 

1.131538 

1.256627 

1.391932 

1.061064 

1.147024 

1.275424 

1.068917 

1.164000 

1.294404 

1.077285 

1.181366 

1.313221 

1 .086387 
1.197909 
1.3313C6 

»* ROW NO. 4 ** 
1.085183 
1.139461 
1.245815 
1.3734L5 

1.091098 

1.148263 

1.260036 

1.388774 

1.097113 

1.158171 

1.273601 

1.402851 

1.103317 

1.169484 

1.266765 

1.415409 

1.109796 

1.182785 

1.304000 

1.116623 
1. 198147 
1.321772 

1.123B24 

1.214631 

1.339654 

1.131554 

1.230742 

1.356976 

** ROW NO. 5 ** 
1.137561 
1.184379 
1.279135 
1.398932 

1.142831 

1.191699 

1.292665 

1.413626 

1.148137 

1.199933 

1.305433 

1.427124 

1.153564 

1.209351 

1.317676 

1.439166 

1.159167 

1.220501 

1.333533 

1.165077 

1.233945 

1.349975 

1.171239 

1.249150 

1.366792 

1.177767 

1.264641 

1.383251 

** MOW NO. 6 ** 
1.190334 
1.230330 
1.313476 
1.424965 

1.194985 
1.236358 
1 .326268 
1.438906 

1.19 9 621 
1.243124 
1.338210 
1.451742 

1.204316 

1.250875 

1.349522 

1.463198 

1.209145 

1.260082 

1.364001 

1.214167 

1.271461 

1.379014 

1.219389 

1.285061 

1.394641 

1.224852 

1.299622 

1.410135 

** ROW NO. 7 ** 
1.243455 
1.277131 
1.348830 
1.451510 

1 .24 7507 
1.282026 
1.360822 
1.464610 

1.251499 
1. 287509 
1.371895 
1.476702 

1.255501 
1.293799 
1.382258 
1.48 7 500 

1.259581 

1.331296 

1.395353 

1.263795 
1. 310662 
1.408875 

1.268150 

1.322435 

1.423200 

1.272657 

1.335736 

1.437623 

*» ROW NO. 8 ** 
1.796B77 
1.324677 
1.385181 
1.478557 

1.300343 
1.328515 
1 .396298 
1.490731 

1.303712 
1.332871 
1.40 6444 
1.501996 

1.307049 
1.337881 
1.415831 
1. 512063 

1.310418 

1.343880 

1.427531 

1.313869 

1.351423 

1.439528 

1.317414 

1.361266 

1.452455 

1.321048 

1.373046 

1.465703 

** ROW NO. 9 ** 
1.350559 

1.353445 

1.356204 

1.358897 

1.361581 

1.364304 

1.367081 

1.369907 




******** 

********* *************** 






*** 



*** 






*** 

ANGLES WITH RESPECT *** 






♦ ** 

TO 

THE l DIRECTION *** 






*** 

OF 

THE ORTHOGONAL 

MESH *** 






* ** 



*** 






******************************** 




** BOTTOM hum NU. 

1 ♦* 








12.563 

13.034 

13.676 


14.484 

15.456 

16.674 

18.303 

20.324 

22.491 

24. 854 

27. 381 


29.279 

29.857 

29.157 

27.917 

26.902 

26. 127 

25.638 

25.47a 


25.650 

25.797 

25.351 

24. 405 

23.233 

21. 841 

20. 220 

18.361 


16.521 





** ROW NO. 2 ** 









1 1 . 4d 4 

11.849 

12.373 


13.051 

13.802 

14.885 

16.247 

17.906 

19. 91 8 

22.070 

24.413 


26.774 

28.046 

28.027 

26.870 

25. 740 

24.871 

24.263 

23.992 


24.067 

24.386 

24.151 

23.341 

22.256 

20.949 

19.410 

17.627 


15.797 





** ROW NO. 3 ** 









10.446 

10.715 

11.132 


11.692 

12.392 

13.229 

14.345 

15.822 

17. 509 

19.442 

21.580 


23.927 

25.842 

26.489 

25.786 

24.561 

2 i. 589 

22.878 

22.496 


22.475 

22.669 

22.846 

22.203 

21. 209 

19.991 

18.539 

16.840 


15.035 





** RUw NU. 4 ** 









9.442 

9.625 

9.944 


10.396 

10.977 

11.662 

12.580 

13.817 

15.265 

16.976 

16.900 


21.046 

23.309 

24.580 

24.519 

23.331 

22.271 

21.477 

20.987 


20.870 

21.251 

21.441 

20.992 

20.093 

18.968 

17.607 

15.998 


14.234 





** MOW NU. 5 ** 









8.465 

8.569 

6.800 


9. 153 

9.625 

10.210 

10.931 

11.950 

13. 180 

14.670 

16.382 


18.316 

20.536 

22.354 

22.903 

22.046 

20. 90 7 

20.040 

19.458 


19.251 

19.547 

19.940 

19.694 

18.908 

17. 881 

16.615 

15.100 


13.394 





** ROW NO. 6 ** 









7.509 

7.542 

7.692 


7.955 

8.326 

8.B01 

9.378 

10.202 

11.235 

12.515 

14.021 


15.746 

17.752 

19. BBS 

20.984 

20.671 

19.486 

18.559 

17.907 


17.613 

17.813 

18.352 

18.307 

17.656 

16. 730 

15.564 

14.14B 


12.514 





** ROW NO. 7 ** 









6.56 9 

6.536 

6.612 


6.792 

7.071 

7.444 

7.9C3 

8.553 

9.410 

10.498 

11.809 


13.333 

15. 125 

17.263 

18.821 

19.046 

17.999 

17.024 

16.326 


15.952 

16.065 

16.686 

16.836 

16.338 

15.517 

14.456 

13.143 


11.594 





** ROW NO. 8 ** 









5.636 

5.547 

5.555 


5.658 

5.852 

6.131 

6.487 

6.984 

7.68 5 

8.601 

9.731 


11.067 

12.654 

14.595 

16.483 

17.178 

16.435 

15.429 

14.701 


14.267 

14.303 

14.953 

15.290 

14.957 

14.24 5 

13.292 

12.087 


10.636 






** ROW NO. 9 ** 



DESCRIPTION OF PROGRAM 


In this section the general characteristics of the program and each of its subroutines 
are given, as well as a dictionary of all the program variables and a complete program 
listing. 

The coding of the program is done entirely in FORTRAN IV. The Lewis library 
routines which are used in the OUPLOT subroutine are described in the appendix. 

Several labeled commons are used throughout the program: 

/INPU/ contains the variables given as input by the user 

/CALC/ contains more calculated constants and the output coordinates and angles of 
the orthogonal mesh 

/CROS/ contains the variables transferred between the MESH and CROSS subroutines 

Several conventions are used in the naming of variables. In any variable name, a 
Z or R usually refers to the Z- and R- coordinates . OM refers to a coordinate or vari- 
able on the orthogonal mesh. BOT and TOP refer to the bottom and top channel walls; 
and LE and TE refer to the leading, or upstream, edge and trailing, or downstream, 
edge of the body, respectively. Slopes have the letters SL in their variable names, and 
SD refers to second derivatives. The letters TEM refer to a temporary storage vari- 
able for some quantity. In the OUPLOT routine a variable containing the letters PLT 
is always an array to be plotted. In the general routines, ROOT, SPLINE, and SPLINT, 
X's and Y’s will always be replaced by Z- and R- coordinates from the calling routine. 

The program is arranged presently to handle a mesh of up to 50 points by 50 points. 

It is being run at Lewis on the IBM 7094/7044 direct coupled system with a 32 768-word 
core (77777^). The total program storage requirement is 47074(g) of which 30000(g) 
is used in storage of variables. The special Lewis plot routines use about 3000(g) words 
of storage. 


Description of Subroutines 

The following sections describe the subroutines which comprise the program. 
Further information on MESH, CROSS, and ROOT can be found in the section DESCRIP- 
TION OF METHOD. 

Subroutine INPUT . - INPUT reads in the input data describing the mesh, the channel 
walls, and the body geometry. It then calculates several constants describing the mesh. 

Subroutine MESH. - MESH is the principal routine in the program. It calculates the 
coordinates of the orthogonal mesh while making use of four other routines - ROOT, 
CROSS, SPLINE, and SPLINT. 
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MESH first estimates the horizontal, or left- to- right, orthogonals by dividing the 
space between channel walls into streamtubes of approximately equal height (see fig. 4). 
These tubes are then intersected by vertical, or bottom-to-top, orthogonal links by 
means of an improved Euler method technique. Each of these links is the average be- 
tween a line which is normal to the bottom orthogonal (of an adjacent orthogonal pair) 
and another which is approximately normal to the top orthogonal (see fig. 5). 

Axial distance between vertical orthogonals is determined by the number of mesh 
lines requested in the following three regions: upstream (to the left ) of the body (ML 
mesh lines from ZBEGIN to ZLE(l)), on the body (MT - ML mesh lines from ZLE(l) to 
ZTE(l)), and downstream of the body (MZ - MT mesh lines from ZTE(l) to ZEND) (see 
fig. 13). MR controls the number of horizontal orthogonals, which is the same in all 
three regions. Note that a vertical orthogonal line which begins within the confines of 
the body may pass across the front or rear edge of the body as the line is generated 
in the radial direction. 

Further details on the operation of MESH are given in the section DESCRIPTION OF 
METHOD. 

Sub rout in e OUTPUT . - OUTPUT prints the coordinates (ZOM,ROM) of the orthogo- 
nal mesh and the angles (AOM) of the horizontal orthogonals with the true horizontal 
plane. These values are printed at each mesh point. 

Subrout ine OUPLOT . - OUPLOT makes two different plots on microfilm for each 
set of input data. The first plot shows the input channel and body geometry. The second 
plot shows the generated orthogonal mesh within the channel region (see fig. 10). 

OUPLOT makes use of six routines from the Lewis in-house library. These six - 
LRMRGN, L RANGE, LRGRID , LRCHSZ, LRLEGN, and LRCURV - are described in the 
appendix . 

OUPLOT first interpolates on the four boundaries to be plotted - bottom and top 
channel walls and front and rear edges of body - to obtain enough points so that these 
boundaries will plot as smooth curves. It then examines the range of values on these 
curves to determine the limits of the two plots. 

OUPLOT then calls the six library routines to accomplish the microfilm plotting. 

A margin is set on the microfilm (LRMRGN), the range of plotted points is given 
(LRANGE), and grid lines and axis labeling are specified (LRGRID). Titles are written 
on the plots by LRLEGN, and finally the desired curves are plotted by calls on LRCURV. 

Subrout ine ROOT . - ROOT is a general routine which locates the point within an 
interval where a function is equal to a given value. In this program the function value is 
RMR (see figs. 7 and 8), and the given value ROOT searches for is RMR = 0. 

ROOT calls another routine (FUNCT) to compute the function of the independent 
variable. In this program CROSS is always called as the function routine, and it calcu- 
lates and returns RMR to ROOT. 
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ROOT begins with a given interval (BNDL to BNDR) in which the intersection 

(causing RMR to equal zero) is known to occur. It then divides the interval in half and 

checks which half contains the intersection. (The check is performed by using the sign 

of RMR. ) The resulting half interval is again divided and checked at the midpoint. By 

o n 

repeating this process 20 times the intersection point is located to within 1/(2 ) of the 

original interval. 

ROOT contains an error message - ROOT HAS FAILED TO CONVERGE IN THE 
GIVEN INTERVAL, which is given if RMR does not approach zero to within the given 
tolerance, TOLER, in 20 iterations. TOLER is set in the MESH routine to 1/2 percent 
of the distance between two adjacent horizontal orthogonals. 

More information on the ROOT routine is given in the section DESCRIPTION OF 
METHOD. 

Subroutine C ROSS. - CROSS is called by ROOT and, in a special case, by MESH. 
When called by MESH, CROSS computes only RCURV (see fig. 8) in the situation where 
the vertical orthogonal link is truly vertical or radial. When called by ROOT, CROSS 
calculates RMR for a given value of the Z- coordinate between BNDL and BNDR. It also 
returns the slope, SL, through ROOT to the MESH routine. 

CROSS uses equations of a cubic spline curve to compute RCURV and SL on the upper 
horizontal orthogonal. The spline curve information is transmitted to CROSS from 
MESH in the /CROS/ common. 

More information on CROSS is given in the section DESCRIPTION OF METHOD. 

Subroutine SPLINE . - SPLINE fits a cubic spline curve to the X- and Y-coordinates 
supplied to it. It solves the tridiagonal matrix equation of reference 2 to obtain coeffi- 
cients for the piecewise cubic polynomial function of the spline curve. Slopes and 
second derivatives at the spline points are returned to the calling subroutine. SPLINE 
uses the end condition that the second derivative at either end point is one-half that at 
the adjacent spline point. 

Subroutine SPLINT . - SPLINT interpolates on a cubic spline curve. It uses the 
same algorithm as SPLINE to fit a curve to the X and Y input points. It then uses the 
interpolation equations of reference 2 to obtain coordinates and derivatives at given in- 
termediate points. 

SPLINT can also be used to extrapolate, although the results are questionable. 

They become more inaccurate the farther the extrapolated point is from the end point 
and the higher the curvature of the spline curve at the end point. When extrapolation is 
done by SPLINT, the message SPLINT USED FOR EXTRAPOLATION is printed, as well 
as coordinates of the spline curve points and extrapolated points. 
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Dictionary of Variables 


The following variables are used throughout the program. The input variables are 
defined in the section Dictionary of input variables. 


A 

AOM 

B 

BNDL 

BNDR 

C 

D 

DELR 

DELZ 

DFX 

DMDZ 

DYDX 

DYDXNT 

D2YDX2 

EOP 

FUNCT 

FX 

FX1 

G 

H 

I 

IPLT 

J 

K 

MAXN 


temporary storage variable in SPLINE and SPLINT 

array of angles made by horizontal orthogonal mesh lines with Z-axis 
at each point of mesh (fig. 9), deg 

temporary storage variable in SPLINE and SPLINT 

left bound on Z- or X- coordinate in MESH and ROOT (figs. 7 and 8) 

right bound on Z- or X- coordinate in MESH and ROOT (figs. 7 and 8) 

temporary storage variable in SPLINE and SPLINT 

temporary storage variable in SPLINE and SPLINT 

incremental distance in R-direction 

incremental distance in Z-direction 

gradient of R with respect to Z along a horizontal orthogonal in ROOT 

expansion distance to size of plot in OU PLOT 

gradient of Y with respect to X in SPLINE 

interpolated gradient of Y with respect to X in SPLINT 

second derivative of Y with respect to X in SPLINE and SPLINT 

end of plot indicator in OUTPLOT 

interval subroutine called in ROOT 

function calculated in ROOT 

function retained in ROOT 

temporary array in SPLINE and SPLINT 

temporary array in SPLINE and SPLINT 

iteration index 

marker indicating which plots have been completed in OUPLOT 
iteration index 
iteration index 

number of spline points or number of interpolated points in SPLINT 
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MEXT 


MINT 

ML 

MLM1 

MLP1 

MR 

MRM1 

MT 

MTM1 

MTP1 

MZ 

MZM1 

N 

N1 

NBOT 

NINT 

NLE 

NTE 

NTOP 

RBOT 

RBPLT 

RBRNG 

RCURV 

REFR 

REFSL 

REFZ 


marker in SPLINT indicating that extrapolation has occurred 

marker in SPLINT indicating that interpolated point coincides with 
given point 

see Input 

ML - 1 

ML + 1 

see Input 

MR - 1 

see Input 

MT - 1 

MT + 1 

see Input 

MZ - 1 

number of X and Y spline points in SPLINE and SPLINT 
N - 1 in SPLINE and SPLINT 
see Input 

number of interpolated X and Y points in SPLINT 

see Input 

see Input 

see Input 

see Input 

array of R- coordinates of points defining bottom boundary of region 
plotted in OUPLOT 

R- coordinate of bottom boundary of points plotted in OUPLOT 

R- coordinate of curved horizontal orthogonal line for a given value 
of Z in CROSS (fig. 8) 

reference R- coordinate of orthogonal mesh point on previously com- 
puted horizontal orthogonal in MESH (figs. 5 and 8) 

reference slope of normal to horizontal orthogonal in MESH (fig. 8) 

reference Z- coordinate of orthogonal mesh point on previously com- 
puted horizontal orthogonal in MESH (figs. 5 and 8) 
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I 


RFPLT 

RLE 

RLINE 

RMR 

RNL 


RNU 


ROM 

RRAD 

RRPLT 

RTE 
RTEM 
RTI, RTI1 
RTOP 
RTPLT 

RTRNG 
SDI, SDI1 
SDTEM 
SL 

SLK 

SLOM 


array of R- coordinates of points defining front, or leading, edge of 
body plotted in OUPLOT 

see Input 

R- coordinate on straight-line vertical orthogonal for a given value 
of Z in CROSS (fig. 8) 

RCURV - RLINE in CROSS (fig. 8) 

R- coordinate of intersection of line normal to previously calculated 
horizontal orthogonal (at REFZ,REFR) with present horizontal 
orthogonal in MESH (fig. 5) 

R- coordinate of intersection of line through previously calculated 
horizontal orthogonal (at REFZ,REFR) with and normal to SLU 
line in MESH (fig. 5) 

array of R- coordinates of mesh points on orthogonal mesh (fig. 9) 

array of R- coordinates of straight vertical lines from bottom bound- 
ary to top boundary of region (fig. 3) 

array of R- coordinates of points defining rear, or trailing, edge of 
body plotted in OUPLOT 

see Input 

temporary storage for R- coordinate arrays 
temporary storage variables in CROSS 
see Input 

array of R- coordinates of points defining top boundary of region 
plotted in OUPLOT 

R- coordinate of top boundary of points plotted in OUPLOT 

temporary storage variables in CROSS 

temporary storage for second derivative arrays 

slope of curve at point of intersection by straight line in CROSS 
(fig- 8) 

temporary storage in SPLINT 

array of slopes at orthogonal mesh points along horizontal orthogonal 
in MESH (fig. 9) 
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SLU 


SXK, SXK1 
SYM, SYN 

TEMP 

TITL1, TITL2 
TITL3, TITL4 

TOLER 

X 

XFIND 

XINT 

XL,XR 

XMED 

Y 

YGIVN 

YINT 

Z 

ZBEGIN 

ZBOT 

ZBOTZ 

ZBPLT 

ZEND 

ZFPLT 

ZLE 


slope of horizontal orthogonal curve where it is intersected by a 
normal from lower orthogonal curve in MESH (fig. 8) 

temporary storage in SPLINT 

data variables in OUPLOT containing identification of symbols to be 
used for plotting 

temporary storage array in MESH 

data arrays containing titles to be plotted in OUPLOT 

tolerance governing calculation of intersection of straight line with 
curve in MESH and ROOT 

array of X- coordinates of spline points in SPLINE and SPLINT 

X- coordinate of intersection of straight line with curve in ROOT 
(fig. 7) 

array of X- coordinates of interpolation points in SPLINT 

left and right X- coordinates of bounds on region in which intersection 
occurs in ROOT 

X- coordinate of midpoint between XL and XR in ROOT 

array of Y- coordinates of spline points in SPLINE and SPLINT 

given value of function, FUNCT, to be matched by root-finding 
process in ROOT 

array of Y- coordinates of interpolation points in SPLINT 

given value of Z-coordinate at which RCURV and RLINE are calcu- 
lated in CROSS (fig. 8) 

see Input 

see Input 

temporary variable in CROSS 

array of Z- coordinates of points defining bottom boundary of region 
plotted in OUPLOT 

see Input 

array of Z-coordinates of points defining front, or leading, edge of 
body plotted in OUPLOT 

see Input 
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ZLRNG Z- coordinate of left boundary of points plotted in OUPLOT 

ZNL Z- coordinate of intersection of line normal to previously calculated 

horizontal orthogonal (at REFZ, REFR) with present horizontal 
orthogonal in MESH (fig. 5) 

ZNU Z-coordinate of intersection of line through previously calculated 

horizontal orthogonal (at REFZ, REFR) with and normal to SLU 
line in MESH (fig. 5) 

ZOM array of Z- coordinates of mesh points on orthogonal mesh (fig. 9) 

ZRPLT array of Z-coordinates of points defining rear, or trailing, edge of 

body plotted in OUPLOT 

ZRRNG Z-coordinate of right boundary of points plotted in OUPLOT 

ZTE see Input 

ZTEM temporary storage for Z-coordinate arrays 

ZTOP see Input 

ZTPLT array of Z-coordinates of points defining top boundary of region 

plotted in OUPLOT 

ZZBOT temporary variable in CROSS 


Program Listing 

A complete listing of the FORTRAN program is given here. If this mesh generation 
program were to be incorporated into a larger program of the reader, the routines 
MESH, ROOT, CROSS, SPLINE, and SPLINT would be used essentially as given here. 
The INPUT and OUTPUT routines would probably be incorporated into input and output 
routines in the reader's program. Finally, the OUPLOT routine would have to be par- 
tially recoded to suit the plotting equipment at the reader’s facility. 


H B FTC MAIN 


10 (.AIL INPUT 
LAI L MESH 
CALL OUTPUT 
0 Al L OUPLOT 

r,u to in 

EN0 
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r> r> 


SlBETC INPUT 


SUBROUT INF INPUT 

c 

0— - 1 NPiiT READS IN INPUT DATA AND CALCULATES SEVERAL CONSTANTS 
C — DESCRIBING THE MESH 
C 

COMMON / l NPU / ML , MT , MZ ,MR»ZBEGIN,ZtND, 

1NBUT , N T OP , ZBOTI50) ,RBOT(50) ,ZTGP(50> .RTOPI50 )• 
?NLE,NTE,7LE< 50) ,RLE(50),ZTE(50> ,RTE<50) 

COMMON /CALC/ MLM1 . ML P 1 • MT M 1 # M T P 1 * MZ Ml * MKM1 * 
l/OM( 50,50 ) # R DM ( 50*50) « AQM ( 5C • 50) 


— READ INPUT DATA 
READ ( 5,1050 ) 
wR I TF ( 6 , 1000 ) 
wR I TEI 6,1050 ) 
W R I TF < 6 , 1 100 ) 
R EAD ( 5,1010 ) 
wRITF(6#10 20 ) 
WRITE! 6,1110 ) 
REaD ( 5,1030 ) 
WR l T E ( 6 • 1 0 40 ) 
wRITE(6,ll 20 ) 
RFAD ( 5,1030 ) 
wRl TF(6,1G40 ) 
WRITF(6,11 30 ) 
RFAD ( 5,10 30 ) 
WRITE! 6,1040 ) 
WRITF(6,11 40 ) 
READ <5,1030 
w R l TF ( 6 , 10 40 ) 
wkfTE(6»l 150 ) 
RFAD < 5,103C ) 
Mk( TE I 6 , 1 040 ) 
WRI TE< 6, 1 160 > 
RFAD (5,1030) 
wRITF<6,10 40 ) 
wRITFl 6,1170 ) 
READ (5,1030 ) 
WK I TF ( 6, 1040 ) 
wRITF(6,118C) 
RFAD (5,1030) 
wRITF(6,10 40 ) 
WRITE(6,119C) 
READ (5,1030) 
wRITF(b,lc4 0 ) 


ML,MT,MZ, MR,NBOT, NT3P*NLE,NTE 
ML# MT , MZ*MR,N BUT, NTQP, NLE, NTE 

ZBEQ I N, Z END 
/BEGIN# ZEND 

< ZBQT ( l ) , 1 = 1 .NBOT) 

(ZBQT ( I )#I=1#NBQT) 

( RBQ T ( I ) . 1 = 1 #NBGT) 

( RBOT ( I ) # 1 = 1, NBG T ) 

( ZTOP< l ) , I =1 ♦ NTQP ) 

( Z T D P ( I ), 1=1, NTQP) 


( RTQP ( [ ) , 1=1 , NTQP) 
( K TOP ( I ) • I = 1 • NT UP ) 

( ZLF( l ) , 1=1, NLE) 

( Z L F ( I ) • I = 1 • NLE ) 

( RLE ( I ) ,1 = 1# NLE ) 

( RLE ( I ) , 1=1 .NLE ) 

( ZTF ( I ) , I = 1 * NTE ) 
<ZTE( I ) ,1=1, NTE) 

( RTE ( I ) ,1=1, NTE) 

( RjE< I ) , 1 = 1, NTE ) 


calculate 

Ml SCtLL A nEQUS CUNSTANTS 

ML Ml 

- 

ML — 1 

ML PI 


ML +1 

MTM1 

= 

MT - 1 

M TP 1 

= 

MT + 1 

M 7 . M 1 

= 

M7-1 

MRMl = 
KFTURN 

MR-1 
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C — FORMAT ST ATEMFNTS 
13UC FORMAT ( 1 Hi I 
1010 FORMAT (16 15) 

1020 FORMAT (41 1X.I6).5X,2I1X,I6>.5X,2(1X,16)) 

1030 FORMAT I8F10.5) 

104C FORMAT <8< 1X.G15.6) ) 

1050 FORMAT (80H 

1 ) 

lion FORMAT ( / / 1 X .3 8FMESH INPUT DATA AND INPUT ARRAY BOUNDS / 5 X , 2HML , 5X . 

12HMT.5X.2HMZ ,5X , 2HMR , 9X , 4H NBOT , 3X , 4HNT0P , 9X. 3HNLE . 4X , 3HNT E ) 

1110 FORMAT (5X.6HZBEGIN.11X.4HZEND) 

1120 FORMAT I / / 1 X . 1 9 H80UNDARY INPUT DAT A/ 5X . 1 C HZ B OT ARRAY) 

1130 FORMAT (5X.1CHRB0T ARRAY) 

1140 FORMAT ( 5 X * 1 GHZ TOP ARRAY) 

1150 FORMAT (5X.1CHRT0P ARRAY) 

1160 FORMAT ( Z Z1X .15HB0QY INPUT DA TA/ 5X , 9HZLE ARRAY) 

117L FORMAT (5X.9FRLE ARRAY) 

1180 FORMAT I5X.9FZTE ARRAY) 

1190 FORMAT (5X.9FRTE ARRAY) 

F NO 


4IBFTC MFSH 

SUBROUTINE MESH 

C 

C — MFSH GENERATES AND CALCULATES THE COORDINATES OF AN ORTHOGONAL 
C — MESH BETWEEN THE INPUT BOUNDARIES 
G 

COMMON / I NP U / ML, MT.MZ. MR. /BEGIN, ZEND, 

1 NBOT, NT OP . Z HOT ( 50) . R BOT ( 5 0 I , Z TOP ( 50 ) ,RT0P(5C ), 

2NL E ,NTE ,ZLE( 50 ) ,RLE( 50) ,ZTE(50 ) ,RTE( 50) 

COMMON /CALC/ MLM 1 . ML PI . M TM 1 , M TP 1 , MZM 1 , MKM 1 , 

1Z0M( 50. 50 ) ,ROM< 50, 50) , AOM( 5C , 5C.) 

COMMON /CROS/ R TEM ( 50 ) , SU TE M ( o( ) , REF Z . REF R , R EE SL 
DIMENSION RR ADI 50.5ul , SLOM< 50 ) , TEMPI 50) 

EXTFRNAL CROSS 

— STORE RRA 0 ON BOTTOM BOUNDARY 

no in 1=1, NBOT 
10 RRAUI 1,1) = RROTI I ) 

— CALCULATE RR AO CN TOP BUUNDARY 

CALL SPLI NT 1 ZTOP ,RTUP »NTUP ,Zb3T , NBOT . RRAU 1 1 .MR ), TEMP) 

— CALCULATE RR AO CN RADIAL LINES FROM BOTTOM TO TCP BOUNDARY 

DO 20 1=1. NBOT 

0 EL R = (RRADII.MR) — RRADI I . 1 ) ) /FLOAT I MRMl ) 

00 2D J =2 . MRMl 

20 RRADI I. J) = RRADII.J-l) +DELR 

C — CALCULATE ZOM ON BOTTOM BOUNDARY 

DEL Z = (ZLEI 1 ) -ZBEG IN ) /FLOAT! MLM1 ) 

/OMIl.l ) = 7 BEGIN 

no 30 I =2, ML 
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30 7 QM (1*1) = ZOM! I-1.1U-UFLZ 

DEL/ = I / TEI 1 ) -ZL £1 1) ) /FLOAT! M T— ML ) 

00 40 I=MLP1.MT 
40 70M1I.1I = ZOM( 1-1, 1 l+DELZ 

OEI Z = (ZEND-ZTE! 1) ) /FLOAT I MZ-MT ) 

00 50 l=MTPl,MZ 
5 t ZOMU.li = ZOM! 1-1,1 H-DELZ 

— CALCULATE ENTRIES IN ZUM AN 0 ROM ROW BY ROW 

00 150 J= 2, MR 

0 

C— CAICUIATE R- COORO I NATE S » SLOPES, AND ANGLES OF PREVIOUS ROW 

CALL SPLINT! ZBOT.RkAD! 1, J-l ). NBOT, ZOM! l.J-l) ,M Z , ROM ( 1 , J-l ) ,SLOM) 
00 60 I =1 , MZ 

60 AOMll.J-1) = ATAN! SLOM! 1) 1*57.295780 

C 

C — CALCULATE RTEM, AND SECOND DERIVATIVES, SDTEM, CF RRAD 
C — ALONG PRESENT ROW 
00 70 1=1. NBOT 
70 K TEM ( I I = RRAD! I, J I 

CALL SPLINE! ZBO T , RTEM , NBOT , TEMP. SDTEM) 

C 

C — MOVE ALONG PRESENT ROW ONE POINT AT A TIME, LOCATING ZOM COORDINATES 
C — OF ORTHOGONAL MESH POINTS ALONG THE ROW 
00 140 1=1. MZ 

REEZ = ZOM! I .J-l ) 

RFER = ROM! I ,J-1 ) 
on 80 K = 1 » NBOT 

IE I7B0T<K).LT.Z0M< l. J-l) ) GO TO 80 

OELR = RRAD! K-l.JI-RRADl K— 1 , J-l I - ! ZOM ( I , J -1 ) — ZBOT ( K-l ) )/ 

II ZBOTI K )-ZBOT(K-l) ) *1 RRAD! K-l . J) -RRAD! K-l , J-l) -RRAD IK, J)* 

2RRADIK, J-l) ) 

GO TO 90 
80 CONTINUE 

90 IF ( ABS < SLUM I I ) ) «L F. G.OuO 1 ) GO TO 120 
C 

C--LOCATF INTERSECTION OF LINE NORMAL TO PREVIOUS ROW WITH PRESENT ROW 
RFFSL = - 1 • / SL 0 M I I ) 

0 El. Z = OELR /RFFSL 
IF (RFESL.LT.O. I GU TO 100 
HNDL = ZOM! I .J-ll 
BNDR = BNOL +2.0*DELZ 
GO TO 110 

ion 8NDR = ZOM! I, J-l) 

BNDL = BNUR+2. 0*DELZ 
1 10 TOLER = DELR/20C. 

CALL ROOT I BN DL . BNDR. 0.. CROSS, TOLER , ZNL ♦ SLU ) 

RNL = REFR+R FFSL*< ZNL-REFZ ) 

GO TO 130 
C 

C — LOCATE INTERSECTION QF LINE NORMAL TO PREVIOUS ROW wITH 
C— PRFSFNT ROW IF NORMAL LINE IS RADIAL 
120 RFFSL =0. 

7 Nl = ZOM! I. J-l) 

CALL CROSS! ZOM! I, J-l I , RNL, SLU) 

C 

C— CALCULATE INTERSECTION OF LINE NORMAL TO PRESENT ROw 
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C — fnOH POINT UN PREVIOUS ROW 

1TL 1 NU = <ZUM! I.J-1HIRUMII ,J-1)-KNL)*SLU+ZNL*SLU**2)/!1.+SLJ**2) 
RNU = RNL-FSLU*! ZNU-ZNL) 

C 

c — calculate final estimate for new zom point by averaging 

140 ZOMII.J) a 1 ZNL+ZNU)/2. 

1 SO CONTINUE 
C 

C— CALCULATE R-COORO I NAT ES , SLOPES, AND ANGLES OF FINAL ROW 

CALL SPLINT! ZBOT , RR AO < 1 , MR ) . NBOT • ZOM ! 1 . MR) ,MZ, ROM! 1, MR) , SLOM) 
DO 160 1=1, MZ 

160 AOM ! I , MR ) = AT AN! SLOM ( I) ) *5 7,2957 8C 

RFTURN 
FNO 


ilKFTC OUTPUT 

SUBROUTINE OUTPUT 

C. 

C— OUTPUT PRINTS THE COORDINATES OF THE ORTHOGONAL MESH 
C 

COMMON / I NPU / ML.MT , MZ . MR , Z BE G l N . Z ENO . 
lNBUT.NT0P,ZB0T(50|,Rb0T< 50) ,Z TOP (501 ,R TOP 150), 
2NLF.NTF ,ZLE ! 50 ) ,RLE1 501 ,ZTE <50 ) .RTF 150) 

COMMON /CALC/ MLM1 , MLP 1 , MTMl , M TP 1 , MZ Ml . MRM1 , 

1 7 OM ( 50 , 50 ) .ROM! 50,50) ,A0M!50,5G) 

C 

C — WKITF ZOM 

WRI TE 16, 1000 ) 
wRITF! 6,1010 ) 

DO 10 J =1 , MR 

IF 1J.FU.1) WRI TE16.1040) J 
IF (J.EO.MR) wR ITE ! 6, 1050 ) J 
IF < J .NE. l.ANO. J.NE .MR ) WR I TE ! 6, 1C 60 ) J 
10 WR I TE! 6.1070 ) ( L OM ( I . J ) . I = 1 , M L ) 

C 

C — WKITF ROM 

wR I TF! 6.1000 ) 
wR I TE! 6 .1020 I 
DO 20 J = 1 , MR 

IF (J.FU.l) wR l TE 1 6 . 1040 ) J 
IF (J.FU.MR) WRITE! 6.1050) J 
IF < J.NE. l.ANO. J.NE. MR) WR I TE 1 6 , 1 0 oO ) J 
20 wRITE(6.l070 ) 1 ROM 1 I.J).I=1«MZ) 

C 

C — WKITF AOM 

w R I T F 1 6 . 1 000 ) 

WRITEI6.1030 ) 

DO 3C J =1 « MR 

IF (J.EU.ll WRITE16.1040) J 
IF IJ.EU.MR) WR ITE! 6. 105C ) J 
IF 1J.NF. l.ANO. J.NE. MR) WR I TE < 6. 1C 60 ) J 
3C wR l TE! 6.1080 I < AOM 1 I , J ) , I = 1 , MZ ) 

RFTURN 



n n no 


c 

C — FORMAT ST ATFMFNTS 
loan FORMAT ( 1 HI ) 

10 1L. FORMAT I / /52X.2 7H** ********************** ***/52X,3H*** ,21 X, 3H***/ 

1 52X.27H*** 7 COORDINATES 3F ***/52X , 2 7H *** ORTHOGONAL MESH 

2 ***/ 5 2X , 3H*** , 21 X, 3H***/52X , 27H*********** ****************) 

10 2 0 FORMAT ( / /52X.2 7H**»************************/52X,3H***,21X, 3H***/ 

1 52X.27H*** R COORDINATES OF *** /52X , 27H *** ORTHOGONAL MESH 

? *** / 1> 2X , 3H*** , 2 1 X , 3H***/52X. 27H*********** ****************) 

10 30 FORMAT I / /50X, 3 2H*************** ********* ********/ 5CX, 3H*** ,26 X, 

13H***/ 5 Cl X , 3 2H*** ANGLES WITH RESPECT * **/5PX , 32H*** TO THE 

2 7 DIRECTION ***/ 50X. 32H*** OF THE ORTHOGONAL MESH ***/5CX,3 
3H***. 26 X, 3H***/ 50X , 32 H** ******************** ********** ) 


10 4C 

FORMAT 

( /IX, 17H** 

BUTTOM 

ROW NO. ,13, 

3H 

10 50 

FORMAT 

< / IX, 14H** 

TOP 

ROW 

NO. , I3.3H 

* * ) 

1360 

FORMAT 

(/IX, 10 H** 

ROW 

NO. 

.I3.3H **» 


10 70 

FORMAT 

( 1X,8G1B.6I 





10 80 

FORMAT 

<1X,8F16*3> 






FNO 







iCHFTC OOPIOT 

SUBROUTINE CUPLOT 
C 

C — OOPLOT PLOTS THE BOUNDARIES OF THE SOLUTION REGION, 

0 — AS v, Ft. L AS THE GENERATED ORTHOGONAL MESH 
C 

COMMON 7 l NPU / ML , MT » MZ » MR .ZBEGIN.ZEND, 

1NBOT . NT OP • ZBOT I 50 I ,RBOT( 50 ) , Z TOP ( 50 ) , R TOP I 5C ) , 

2NL F «NT F ,2 L E I 50 I ,RLE I 50 ) , Z TE I 50 ) , R TE ( 50 I 
COMMON /CALC/ MLM1 . ML P 1 , MTM1 , MTP1 . MZM1 ,MRM1 , 

1/0M( 50, 50 > .ROM! 50, 501 ,AOM( 5C, 5G> 

0 IMFNSION ZBPLTI 13C ) , RBPL T 1 10 0 ) ♦ Z TPLT ( 1 03 I ,R TPLT ( IOC), 

1/ FPLT < luO ) , RFPL T ( 100 ) , ZR PL T 1 1 0(. I , RRPLT ( ID 0 I , ZTEMI1O0 ) ,RTEM( ISC ) , 
2SQTFM< 1G0 ) ,T ITL 1(5 ) .TITL2I3) .TITL3I2 I .TITL4I 21 
DATA TITI1/30H BOUNDARIES OF PLOT REGION / 

OAT A TITL2/18H ORTHOGONAL MESH / 

DATA TITL3/12HZ DIRECTION/ 

DATA TITL4/12HR DIRECTION/ 

OATA SYM/IHX/ 

DATA SYN/1H0/ 

— OBTAIN PI OT POINTS ON LOWER BOUNDARY 
DEL / = (ZBOT (NBOT ) —ZBOT ( 1) ) / 99 • 

/ BPL T ( l I = ZB0TI1I 
ZBPLT(IOO) = ZBOTINBOT) 

DO 10 1=2,99 

10 ZBPLTIII = ZBPLTI l-ll PUFLZ 

CALL SPLINT I ZBO T , RBOT , NBOT • ZB PL T, 103 , KBPL T , S OT EM ) 

— OBTAIN PLOT POINTS ON UppER BoUNOarY 
3ELZ = ( 7 T OP IN TOP ) — Z T OP (1)1/99. 

ZT PI Till = Z TOP 11) 

ZTPLTI100I = ZTOPINTOPI 
DO 20 1=2.99 
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20 7TPLTII) = ZTPLTI l-ll+OELZ 

CALL SPL I NT < ZTQP.RTUP.NT OP . ZTPL T . 100 . R TPL T , S DT EM I 

— DETAIN PLOT POINTS UP LEADING EDGE OF BODY 
DELR = IRLEI NLE I-RLEI 111/99. 

KFPLT1 II = R I F I 1 1 
RFPLTI 100 ) = RLE(NLE) 

DO 30 J =2 » 99 

30 rFPITIJI = RfPLTI J-1I+DELR 

CALL SPLINT! RLE ,ZL E . N LE . R FPL T . 100 . ZF PL T , SDTE MI 

— OBTAIN PLOT POINTS UP TRAILING EDGE UF BODY 
DEL R = (RTEINTE l-RTFI 1) 1/99. 

RRPLTI 1 > = RTEI II 
RKpLT I 100 I = RTEINTE) 

DO AO .1=2 « 9 9 

AO RRPLTIJ) = RRPLTI J-1M-DELR 

CALL SPLINT! RTE ,Z T E . NTE . RRPL T , 100, ZRPLT, SDTE MI 

— CALCUI ATE THF RANGE OF THE PLOT 

Zl.kNG = AMIN1I ZBOTI 1 I tZTOPI 1) 1 
7 RRNG = AMAXK ZBOTINBOTI .ZTOPI NTUP) I 
DEI Z = ZRRNG-Z 1 RNG 
ZLRNG = 7 L RNG— 0 .05*DELZ 
7 RRNG = 7RRNG+O.Q5*OELZ 
R RRNG = R HOT I 1 I 
DO 50 I =2 » NB OT 

51. RBRNG = AMINKRBRNG.RBOT! II I 
KTRNG = R TOP 1 1 ) 

DO 60 I =2 • N TOP 

60 RTRNG = AM AX 1 1 RTRNG * R TOP I I I I 
DELR = KTRNG-RBRNG 
kBRNG = RBRNG— 0.0 5# DELR 
RTRNG = RTRNG+0.05*DELR 

— CHOOSE MAXIMUM RANGE. AND EXPAND RANdE IN OTHER DIRECTION 

DMD? = 1. 1*ABS( DELZ-DELR) /2. 

IF I DELR. GT.OELZ I GO TO 70 
RTRNG = RTKNG+0M02 
RBRNG = KBRNG— UMD2 
GO TO BO 

7( 7 RkNG = 7RRN6+DMD2 

7.1 RNG = ZLRNG-DMD2 
f. 

C — PLOT BOUNDARIES OF REGION AND PL J T ORTHOGONAL MESH 
BO CALL I KMRGNI 1.0,1.0.1.0,1.01 

CALI I.RANGEI ZL R NG. Z RRNG, RBR NG , RTRNG I 
CALL I RGR ID! -1. -1,1.0. 1.0 I 
I PI T = 1 

CALL LRLFGNI TITL1. 2 8,0.3.0,C. 5,0.0 
9D IF IIPLT.FD.2I CALL LRLEGNI T I TL2 . 1 6 , 0 , 3 . 6 ,0 . 5, 0. C I 
CAI I. I RCHS7 I 21 

CALL LRLEGNI T I T L3 , 1 2 , 0 . 7. 5 . C . 8 , C . 0 I 
CALL LRLEGNI Tl TLA, 12,1 ,C. 8, 7. 5, 0.0 I 
0 At L LRCHSZ I 31 

CALL I Kf.URV! ZBPLT « RBPLT. IOC »2. SYM.L.C I 
CALI LRCUKV! ZTPLT.RTPLT « luG , 2 , SYM.C.O) 

CALL I RCURVI ZFPLT.RFPLT. 10L .2, SYM.C.O) 

CALL I RCURV I ZRPLT .RRPLT. IOC, 2. SYM.C.O ) 
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IF IIPLT.E0.2) GO TO ICG 

CAl L LRCURVI 7BGT.RBOT,NBOT,4,SYM,0.0) 

CALL I RCURVI 7TOP , RTQP , NTGP , 4 , S YM ,C ,Q ) 

CAl I LRCURVI ZLE. RLE, NLF, 3, SYN. 0.0) 

CALL LRCURVI 7T E ,R T E . NTE . 3 , S YN . 1. 0 ) 

I PL T = 2 
GO TO 90 

C 

C— PLOT VERTICAL MESH LIMES 
100 on 120 1=1, M7 

7TEMI 1 l = 70 Ml 1,11 
RTFMI U = RO Ml 1,11 
DO 110 J= 2 « MR 
7TEMI .J l = 20MI I , J I 
11L RTFMI J) = ROM! I , J ) 

1 2t. CALL I RCURVI ZTEM.RTEM.MR ,2 .SYM.C.C) 

C 

C — PLOT HORIZONTAL MESH LIMES 
FOP = 0,0 
DO 130 J= 2 . MRM1 
IF IJ.EU.MRM1) EOP= 1,0 

130 CALL LRCURVI 70MI1,J) , ROM I 1 , J I , MZ,2,SYM,E0P1 
RETURN 
FMD 


ilBFTC ROOT 


SUBROUT INE ROOT lBNDL,BNOR,YGI VM, FUMCT , TOLER, XF INU, OF XI 

C. 

f. — ROOT FINOS A ROOT FOR I FUNCT-YGI VN) IN THE INTERVAL (BNDL.BNOR) 
C 

XL = BNDL 
XR = BN DR 

CALL FIJNCTI BNDL ,FX1 ,OFX) 

DO 2D 1=1,20 
XMID = I XL + XR) / 2, 

CAl L FU NO T ( X M I D * F X . OF X I 

IF I ( FX1-YGI VNl*< FX-YGI VN) .GT.C. I GO TO ID 
XR = XMID 
GO TO 20 
10 XL = XMID 
F X 1 = FX 


2D CONTINUE 

XFIND = XMID 

IF I ABSlYGIVN-FXI.LT. TOLER) RETURN 
wRITE(6,1 OGC ) BNDL.BNDR.YGIVN 
STOP 

1000 FURMATI // //4X.49HRU0T HAS FAILED TO CONVERGE IN THE GIVEN INTERVAL 
1 7 R-X , 6HBND L = ,G 1 4. 6, 10X . 6HBMDR = , G1 4. 6 , 10 X , 7H YG i VN =,G14.6) 

F NO 
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SIBFTC CkOSS 


SUBROUTINE C ROS S < Z ,RMR ,SL ) 

f. 

C — r.RUSS CALCULATES R AS A FUNCTION OF Z ALONG A CURVE AND A 
C — STRAIGHT LINE WhlCH INTERSECTS IT. AND THEN FOR A GIVEN VALUE OF Z 
C — CALCULATES THE DIFFERENCE BETWEEN THE VALUES OF R ON THE 
c — STRAIGHT linf and the curve 
c 

COMMON /I N P U / ML . MT , MZ . MR , Z BEG IN . Z END . 

INBOT.NTOP . ZBOT <501 ,RB0T(5D) »ZT0P(5G) ,RT0P(5Q I, 

2NLE ♦ NT E *Z L E ( 5C » .RL E < 50 > • Z TE < 50 1 . R TE < 50 » 

COMMON /CROSZ R TEM < 50 ) , SDTE M< 5C ) , REF Z . k EF R , R EF SL 
C 

C— LOCATE POSITION OF Z IN RRAD (OR RTEM ) ARRAY 
00 10 I =2 • NBOT 
IF (Z.LF.ZBOTUn GO TO 20 
10 CONTINUE 
C 

C — COMPUTE R-CDORD IN ATE (RCURV) AND SLOPE ( SL » ON THE CURVE 
C — FOR THE GIVFN VALUE OF Z 

20 0 EL Z = ZBOT( II-ZBOTI 1-1) 

7. ROT 7 = /BOTH >-Z 
ZZBOT = Z - Z B OT < I-ll 
SOI = SOT EM ( I t 
SDH = SUTEM ( I — 1 ) 

RT I = RTEM< I I/DELZ 
kT I 1 =■ RT EM ( 1-1 ) ZOELZ 

k C UR V = S0ll*ZB0TZ**3/6./DELZ *■ SDI*ZZB0T**3/6./DELZ *■ 

1< RT l-SDI*DELZZ6. )*ZZBOT ♦ < RT I 1- S D 1 1 *OELZ /6 . » * ZBOT Z 
SL = -SDI1*ZB0TZ**2/2.ZDELZ + SOI *ZZBOT**272 .7DELZ + 

IRTI-kTIl - ( SDI-SD ll»*DELZ/6. 

IF 4RFFSL.EQ.0. ) GO TO 30 
C 

C — COMPUTE R-COORDINATE (RLINE) ON STRAIGHT LINE 
C — FOR GIVEN VALUE OF Z 

RL INF = RffH *REFSL*( Z-REFZ ) 

C 

C — COMPUTE DIFFERENCE IN R COORDINATES 
R MR = KCURV-RLINE 
kFTURN 
C 

C — SPECIAL CASF FOR RADIAL STRAIGHT LINE 
30 RMR = RCURV 
RETURN 
FND 


MBFTC SPLINE 

SUBROUT INF S PL I NE < X . V , N. 0 VOX , D 2 YDX2 ) 

C 

C — SPLINE FITS A SPLINE CURVE TO X ANU Y 

C--ANO CALCUl ATES FIRST AND SECOND DERIVATIVES AT THE SPLINE POINTS 
C — END POINT SECUND DERIVATIVES LNE HALF THOSE AT ADJACENT POINTS 



f) IMF NS ION X I N) * YIN) * OYDX ( N ) , D 2 YDX2 ( N > 

u i mens i on (;uoni,Nd(0) 

s < n = -0.5 

H { 1 I = 0 . 

*\i L = N- 1 

IE C Ml • LT • 2 I G( 11=1. f 
IE <N1.LT .21 GO TO ?r 
00 1^ I =2 * Nl 
A = (Xi \ >-X< i~l) ) / 6. 
rt= UUf 1 )-X < I ) >/6. 

C = ?.*( A^I-A*C,U-I1 

0 = (YU+ii-vini/iximi-xim-iYin-yu-iii/ixm-xu-ii) 

G ( l ) = »/C 

If rti n= ( f)-A*H< I- II ) /c 
? r D2Y0X2INI = HINH/I2.0+GIN1)) 

1)0 3D I =2 * N 
X = N + l- I 

*0 .)2YDX?<Kl = H<M-G<Kj*D2YDX2<K+li 

OYDX< U= i XI 1)-X< 2)) /6.*< 2.*D2Y0X2 ( 1 ) +D2YDX2 (2 ) ) f ( Y(2) -Y ( 1) ) / ( X( 2) 
l-Kun 
NO AO l=2.N 

4C 0 Y OX <11= (X<n-X<I-m/6.*( 2.*02YGX2( I ) + D2 YGX2I t-1 ) IMYU |-Y( 1-1) ) 
1/ < XI f >-X< 1-1 ) ) 
nETIJRN 
END 


% IrtETC SPL 1 NT 

SUBRoUT IN E S PL I NT < X • Y ♦ N * X I NT , N I NT , Y I N T , DY 0 XN T) 

C 

C— SKINT FITS A SPLINE CURVE TO X ANU Yt AND THEN CALCULATES 
C. — INTER POL AT FO Y VALUES ANU DERIVATIVES AT GIVEN X POINTS 
0~ ENU POINT SECOND DERIVATIVES ONE HALF THOSE AT ADJACENT POINTS 
C 

DIMENSION X I M • Y(N) ♦ X I NT I N l N T I tYINT(NINT) * DY DXNT ( N i NT ) 

U l MENS I ON G ( inO)«HC10C;ltD2YDX2ilC0) 

IE (NINT.LE.C) RETURN 

C 

C COMPUTE SPLINE CURVE 
G m = -0 . s 
h ( i i = o. 

Nl= N— l 

IF <N1 . LT . 2 I G ( 11*1.0 
[ f <Nl.LT.2I GO TO 20 
00 10 I =2 * N 1 
A — u< I )— X I I -II ) / 6. 
rt= IXUfli-XlI) )/6. 

C= 2.*I A+ft)-A*G< 1-11 

0= < y< i+d-yi m/ixu + n-xi m-< yih-yi i-in /< xi n-xi i-in 

iiin= b/c 

in H< i 1= (0-A*H(I-l) 1 /c 
2D D2YOX24N) = hi N 1 1 /( 2. O +G { N 1 } I 
on 30 l =2 • N 
K- N +■ 1 - I 
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30 02YDX2IK1= HIKl-G(K)*Q2YDX2<KFl) 

f. 

f. I NTFRPOLATE ON SPLINE CURVE TO GET YINT 
MINT = ft 
ME XT = 0 
DO 120 1= 1 . N INT 

K = 2 

IE (XINTIII-Xmi 50.40.60 
40 Y INTI H= YI 1 ) 

MINT = 1 
GO TO 100 
50 MF XT = 2 
GO TO 100 

60 IF (XINTI l 1-XlK) 1 100,70,80 

70 Y INTI 11= YIK » 

MINT = 1 
GO TO 100 
80 K = KFl 

IF (K-Nl 60.60.90 
90 <= N 

MFXT = 2 

100 SLK = XIK1-XIK-1I 
S X K = X I K> -X I NT I II 
SXKl = X INT I I >-X IK-1 » 

IF IMINT.EO.il GU TO 11C 

Y INTI I 1 = D2YDX2I K-l ) *SXK**3/6. /SLK+D2 Y0X2 ( K 1 *S XK. 1 * *3 /6 . / S L K 
1#-I Y l Kl / SLK-D2YGX2! K»*SLK/6. 1* SXK1* ! Yl K- 1 1 / SL K-02YDX2 IK-1 > *SLK/fc. 1 
2*SXK 

1 1 C MINT = 0 

DYDXNTI 11= — 02 Y 0X2 I X- 1 1 *SXK**2 /2./SLKFD2YDX2 !K)*SXKl*«‘,i/2./SLK 
1*1 YI K 1-Y< K-l 11/ SL<- I 02Y0X2 I K 1 -D2YOX2 I K-l 1 ) *SLK/6. 

IE IMFXT.E0.2l WRITE! 6,1000 1 X INT I I 1 . Y I NT I I ) 

IF IMFXT.EO.?) ME XT = 1 
1 2u CONTINUF 

MAXN = MAXO(N.NINT) 

IF <MFXT.GT.G> WRi TE I 6. 101C ) N.NINT, I XI l ) , Y 1 1 1 . X I NT ( I 1 . Y INTI I ) , 

1 DYDXNTI I » , I =1, MAXN I 
R ETURN 

lOOu FORMAT ( 56 H SPLINT USEO FOR EXTRAPOLATION. EXTRAPOLATED X ANO Y = 
1. 2G14.6 1 

IDIO FORMAT! 2X . 24FN0MHER OF PUINTS GIVEN = , I 3 . 4 X . 31 HNUMB ER OF POINTS IN 
ITER POL A TF D = . I 3 / 10X . 1HX , 1 9X , H Y . 14X. 1 4HX l NT ER POL A T E D , 6 X , 1 4 HY INTE 
2RP0I AT FO. 5 X . 17HDYDX INTERPOLATED/! 5E2C. 8) 1 
END 


CONCLUDING REMARKS 

A FORTRAN IV program has been presented which computes and plots coordinates 
for a two-dimensional orthogonal mesh between the walls of a channel. All routines in 
the program are in easily transferable FORTRAN code. A plot routine, however, uses 
in-house subroutines and equipment and would have to be recoded elsewhere. 


\ 
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This report has presented a description of the method used, a description of the 
program itself, and a complete listing. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, January 28, 1972, 

764-74. 


38 



APPENDIX - LEWIS LIBRARY SUBROUTINES 


The following descriptions of the NASA Library subroutines - LRMRGN, LRANGE, 
LRGRID, LRLEGN, LRCHSZ, and LRCURV - are copied from the NASA Lewis Program- 
mers Manual. These routines are part of a microfilm plotting package called CINEMATIC. 


Subroutine LRMRGN 

Purpose. - LRMRGN is used to change the width of plot margins. 

Usage. - CALL LRMRGN (XLEFT, X RIGHT, YBOTM, YTOP) where XLEFT (float- 
ing point) is the left margin width in absolute positioning units, XRIGHT (floating point) 
is the right margin width in absolute positioning units, YBOTM (floating point) is the 
lower margin width in absolute positioning units, and YTOP (floating point) is the upper 
margin width in absolute positioning units. 

Method , - CINEMATIC sets margins around the plotting area as follows: LEFT and 
BOTTOM, 1. 0 absolute positioning units; RIGHT and TOP, 0. 4 absolute positioning 
units. A call to LRMRGN before LRCURV will change the width of the margins. A 
frame of film contains 10 absolute positioning units in the horizontal and in the vertical 
directions. Thus, a margin of 1.0 absolute positioning units is 1/10 of a frame wide. 

Restrictions . - LRMRGN must be called before the plotting routine it applies to. 

The settings of LRMRGN remain in effect until changed by another call to LRMRGN. 


Subroutine LRANGE 

Purpose. - LRANGE is used to set the range of (X,Y) curve points. 

Usage. - CALL LRANGE (XLEFT, XRIGHT, YBOTM, YTOP), where XLEFT is the 
left end point of a plot in the user’s units, XRIGHT is the right end point of a plot in the 
user's units, YBOTM is the lower end point of a plot in the user's units, and YTOP is 
the upper end point of a plot in the user's units. 

Method . - The curve-plotting subroutine LRCURV searches the (X,Y) coordinates 
for maximums and minimums and scales the rest of the user’s points to fit between 
them. A call to LRANGE before LRCURV suppresses the search. The settings of 
LRANGE remain in effect for all successive plots until changed by another call to 
LRANGE. XLEFT = XRIGHT = 0. 0 can be used to remove the LRANGE X-settings with- 
out providing new ones. YBOTM = YTOP = 0.0 does the same for Y-settings. 
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Restrictions . - LRANGE must be called before the curve plotting routine it applies 
to. The settings of LRANGE remain in effect until changed by the user. 


Subroutine LRGRID 

Purpose. - LRGRID is used to specify grid line changes. 

Usage. - CALL LRGRID (IXCODE, IYCODE, DX, DY). IXCODE (fixed point) is a 
switch used as follows: It applies to vertical grid lines. 

IXCODE = 0 means return to using CINEMATIC'S built-in grid format (11 grid 
lines). 

IXCODE = ±1 means DX specifies how many grid lines; IXCODE = -1 suppresses 
grid labels. 

IXCODE = ±2 means DX specifies grid intervals; IXCODE - -2 suppresses grid 
labels. 

IXCODE = ±3 means DX specifies how many "tick marks" instead of grid 
lines; IXCODE = - 3 suppresses grid labels. 

IXCODE = ±4 means DX specifies interval between "tick marks"; IXCODE = -4 
suppresses grid labels. 

DX (floating point) specifies grid line or "tick mark" frequency or intervals, depending 
on how IXCODE is set. IYCODE (fixed point) is the same as IXCODE, but it applies to 
horizontal grid lines. DY (floating point) is the same as DX but for horizontal grid lines. 

Method . - CINEMATIC puts 11 horizontal and 11 vertical grid lines on every plot, 
unless LRGRID is called. When a grid line frequency is specified, CINEMATIC sets the 
interval between the specified number of grid lines to be equal to ZxlO n , where 
Z = 1. 0, 2.0, 2.5, or 5. 0 and n depends on the magnitude of the user's data. To get 
these intervals, CINEMATIC will adjust the end points of the plot, if necessary. This 
adjustment occurs only when a grid line frequency is specified. To avoid any adjust- 
ments, specify grid intervals. 

Restriction . - LRGRID must be called before the plotting routine it applies to. The 
settings of LRGRID remain in effect until changed by another call to LRGRID. 


Subroutine LRLEGN 

Purpose, - LRLEGN is used to print a legend anywhere on a plot. 

Usage. - CALL LRLEGN (CHARS, N, IORIEN, X, Y, EOP). CHARS is an array 
of characters to be printed. CHARS must be dimensioned large enough for the number 
of characters desired (four characters per word on the 360; six characters per word 
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on the 7094/7044 direct couple). N (fixed point) is the number of characters to be 
printed. IORIEN (fixed point) is a switch: 

IORIEN = 0 causes horizontal printing. 

IORIEN = 1 causes vertical printing. 

X (floating point) is the X- coordinate of the starting point in absolute positioning units. 
Y (floating point) is theY- coordinate of the starting point in absolute positioning units. 
EOP (floating point) is a switch: 

EOP = 0 indicates the current plot is not yet complete. 

EOP = 1 indicates the current plot is complete. No more calls to plotting or 
printing subroutines for this plot will occur. 

Method. - The user expresses the (X,Y) starting point of a line of printing in ab- 
solute positioning units. Absolute positioning units range from 0. 0 to 10. 0 in both the 
X- and Y-directions for one frame of film. Absolute positioning units give the user a 
coordinate system to specify (X,Y) points independent of points on a curve. 

LRLEGN prints medium-size characters. The user may also get other character 
sizes, italics, lower case, and special symbols. 


Subroutine LRCHSZ 

Purpose. - LRCHSZ is used to change the size of printed characters. 

Usage. - CALL LRCHSZ (ISIZE) where ISIZE (fixed point) gives the size: 

ISIZE = 0 means let CINEMATIC resume selecting the size. 

ISIZE = 1 means miniature characters. 

ISIZE = 2 means small characters. 

ISIZE = 3 means medium characters. 

ISIZE = 4 means large characters. 

Method. - LRCHSZ changes the character size for all character printing that fol- 
lows. The specified size remains in effect until changed by another call to LRCHSZ. 
Large: 43 characters per line, 22 lines per frame. 

Medium: 64 characters per line, 32 lines per frame. 

Small: 86 characters per line, 43 lines per frame. 

Miniature: 128 characters per line, 64 lines per frame. 

Restriction . - LRCHSZ must be called before the printing subroutine it applies 
to. 


Is. 
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Subroutine LRCURV 


Purpose . - LRCURV is used to plot one curve of a multiple- curve plot. 

Usage. - CALL LRCURV (X,Y,N, I TYPE, SYMBOL, EOP). X (floating point) is 
an array of X- coordinates for the curve. Y (floating point) is an array of Y- coordinates 
for the curve. N (fixed point) is the number of (X, Y) points to be plotted. ITYPE is a 
switch that indicates the type of plot desired: 

ITYPE = 1 specifies a dot plot. Each (X, Y) point is represented by a dot. 

ITYPE = 2 specifies a vector plot . Successive (X,Y) points are joined by 
straight lines. 

ITYPE = 3 specifies a symbol plot . Each (X, Y) point is represented by a sym- 
bol. The FORTRAN character in SYMBOL specifies the symbol used. 

ITYPE = 4 specifies a special symbol plot. Each (X, Y) point is represented by 
a special symbol taken from the SPECIAL CHARACTER TABLE. The 
special symbol used is the one corresponding to the FORTRAN character 
in SYMBOL. 

ITYPE = 5 is the same as ITYPE = 3, except that a smaller size symbol is 
used. 

ITYPE = 6 is the same as ITYPE = 4, except that a smaller size symbol is 
used. 

SYMBOL specifies the plotting symbol when ITYPE = 3 or 4. When ITYPE = 1 or 2, 
SYMBOL must appear in the call list but is not used by LRCURV. EOP is a switch that 
indicates when the last subroutine call for a given plot is being made: 

EOP =0.0 means the current plot is not yet complete. More subroutine calls 
for this plot will follow. 

EOP =1.0 means the current plot is complete. No more printing or plotting 
subroutines will be called for this plot. 

Method , - LCURVE provides greater flexibility in drawing curves. LRCURV is 
useful for the plotting situation in which all (X,Y) points for a plot are not in the com- 
puter memory at the same time. Several calls to LRCURV may be made for the same 
plot. 

The X and Y arrays are in whatever units the user is working with. LRCURV scales 
his data range to fit the size of the plot on film. The user should call LRANGE before 
LRCURV to supply the range of his data points to CINEMATIC. If the user does not 
call LRANGE, LRCURV will take the user’s data range from the first call to LRCURV 
for any given plot. 

LRCURV does not destroy the contents of X, Y, N, ITYPE, SYMBOL, or EOP dur- 
ing plotting. 
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